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Abstract 

We give an explanation for the onset of wind-induced flutter in a flag. Our theory accounts 
for the various physical mechanisms at work: the finite length and the small but finite bending 
stiffness of the flag, the unsteadiness of the flow, the added mass effect and vortex shedding from 
the trailing edge. Our analysis allows us to predict a critical speed for the onset of flapping as well 
as the frequency of flapping. We find that in a particular limit corresponding to a low density fluid 
flowing over a soft high density flag, the flapping instability is akin to a resonance between the 
mode of oscillation of a rigid pivoted airfoil in a flow and a hinged-free elastic filament vibrating 
in its lowest mode. 
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The flutter of a flag in a gentle breeze, or the flapping of a sail in a rough wind are 
commonplace and familiar observations of a rich class of problems involving the interaction 
of fluids and structures, of wide interest and importance in science and engineering 
Folklore attributes this flapping instability to some combination of (i) the Benard- von 
Karman vortex street that is shed from the trailing edge of the flag, and (ii) the flapping 
instability to the now classical Kelvin-Helmholtz problem of the growth of perturbations at 
an interface between two inviscid fluids of infinite extent moving with different velocities . 
However a moment's reflection makes one realize that neither of these is strictly correct. The 
frequency of vortex shedding from a thin flag (with an audible acoustic signature) is much 
higher than that of the observed flapping, while the initial differential velocity profile across 
the interface to generate the instability, the finite flexibility and length of the flag make 
it qualitatively different from the Kelvin-Helmholtz problem. Following the advent of high 
speed flight, these questions were revisited in the context of aerodynamically induced wing 
flutter by Theodorsen jj| . While this important advance made it possible to predict the onset 
of flutter for rigid plates, these analyses are not directly applicable to the case of a spatially 
extended elastic system such as a flapping flag. Recently, experiments on an elastic filament 
flapping in a flowing soap film (|, and of paper sheets flapping in a breeze |f| have been used 
to further elucidate aspects of the phenomena such as the inherent bistability of the flapping 
and stationary states, and a characterization of the transition curve. In addition, numerical 
solutions of the inviscid hydrodynamic (Euler) equations using an integral equation approach 



and of the viscous (Navier-Stokes) equations £| have shown that it is possible to simulate 
the flapping instability. However, the physical mechanisms underlying the instability remain 
elusive. In this paper, we aim to remedy this using the seminal ideas of Theodorsen 

We will start by considering the dynamics of an inextensible one-dimensional elastic 
filament of length L and diameter d and made of a material of density p s and Young's 
modulus E embedded in a two dimensional parallel flow of an ambient fluid with a density 
p and kinematic viscosity u, shown schematically in Fig. HJ 1 . We assume that the leading 
edge of the naturally straight filament is clamped at the origin with its tangent along the 
x axis, and that far from the filament, the fluid velocity U = £/x. Then the transverse 



1 Our analysis also carries over to the case of an elastic sheet a 3-dimensional parallel flow with no variations 
in the direction perpendicular to the main flow. 
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FIG. 1: Schematic representation of the system. An elastic filament of length L, clamped at the 
origin is embedded in a 2-dimensional flow of an inviscid fluid with velocity U in the x direction. 
Its lateral position is denoted by Y(x,t). 

position of the filament Y(x,t) satisfies the equation of motion 1.0]: 

mY tt = -BY XXXX + IAP. (1) 

Here, and elsewhere Af, = dA/db, m = p s ird 2 /4: is the mass per unit length of the filament, 
B = it Ed 4 / 64 its flexural rigidity, I is the thickness of the fluid film 2 and AP the pressure 
difference across the filament due to fluid flow. In deriving (JTJ) we have assumed that the slope 
of the filament is small so that we can neglect the effect of any geometrical nonlinearities; 
these become important in determining the detailed evolution of the instability but are not 
relevant in understanding the onset of flutter. For the case when the leading edge of the 
lag is clamped and the trailing edge is free, the boundary conditions associated with are 



lag i 

m 



Y(t,0)=0, Y x (t,0) = 0, 
Y xx (t,L)=0, Y xxx (t, L) = 0. (2) 

To close the system (1 1121) we must evaluate the fluid pressure AP by solving the equations 
of motion for the fluid in the presence of the moving filament. We will assume that the flow 
is incompressible, inviscid and irrotational. The omission of viscous effects is justified if the 
shear stress pJ induced by the Blasius boundary layer [l2| is small compared to the 
fluid pressure pU 2 far away from the filament or equivalently if the characteristic Reynolds 
number Re = UL/u ^> 1. In typical experiments, since Re ~ 10 5 , this condition is easily 



2 In the experiments with filaments in soap films -Jo], I ^ d. For a sheet I = 1, m is a mass per unit area 
and B is now the bending stiffness per unit length. 
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met. Then we may describe the unsteady fluid flow as a superposition of a non-circulatory 
flow and a circulatory flow associated with vortex shedding, following the pioneering work 
of Theodorsen j^J. This allows us to respect Kelvin's theorem preserving the total vorticity 
of the inviscid system (which is always zero) by considering a vortex sheet in the fluid and 
an image sheet of opposite strength that is in the filament. Both flows may be described 
by a velocity potential <fi which itself may be decomposed into a non-circulatory potential 
(j) nc and a circulatory potential 7 with <ft = <p nc + <p 1 . Then <ft satisfies the Laplace equation 
V n 20 = characterizing the two-dimensional fluid velocity field (u, v) = (<f) x ,<j)y). 

For small deflections of the filament, the transverse velocity of the fluid v varies slowly 
along the filament. Then we may use a classical result from airfoil theory for an airfoil 
moving with a velocity v — Y t + UY X to deduce the non-circulatory velocity potential along 
the filament as 



nc = y/x (L-x) [Y t + UY X ] , (3) 

To determine the jump in pressure due to the non-circulatory flow we use the linearized 
Bernoulli relation so that 

P nc = -2p(d t (f) nc + Ud x (f) nc ) (4) 
pU(2x - L) 



/ (Yt + UY X ) + y / x(L- x)pY a . 

y x(L — x) 

Here we note that the fluid added-mass effect 3 is characterized by the term proportional 
to Y tt , and we have neglected terms of order 0(Y xt ) and higher associated with very slow 
changes in the slope of the filament. 

Kelvin's theorem demands that vorticity is conserved in an inviscid flow of given topology. 
Thus, the circulatory flow associated with vortex shedding from the trailing edge requires a 
vorticity distribution in the wake of the airfoil and a (bound) vorticity distribution in the 
airfoil to conserve the total vorticity. If a point vortex shed from the trailing edge of the 
filament with strength — T has a position ^ (1 + X ), X Q ^ 1, we must add a point vortex 
of strength T in the interior of the sheet at |f ^1 + -^J. This leads to a circulatory velocity 
potential along the filament 

T / yJx(L - X) V / X 2 - 1 

(Pr — arctan 1 



27T ^ |(1 + x ) - xx 



When the filament moves, fluid must also be displaced and the sheet behaves as if it had more inertia 



where Xq = x " + ^/ Xo characterizes the non-dimensional center of vorticity which is at (1 + 
Xq)/2. Therefore for a distribution of vortices of strength 7 defined by T = "f^dx , the 
circulatory velocity potential is 

il r ( y^y^ix 



07 = -77-^7 / arctan — lTTT X 7dXo ' 

2-n2 J 1 y 1 (1 + x ) - xx J 

To calculate the pressure difference due to the circulatory flow, we assume that the shed 
vorticity moves with the flow velocity U in the flow so that <9 4 7 = jUd x ^^ 4 . Then, we 
may write j^j: 

p., = r 2 * + ^_-y (6) 

2ix^/x{L — x) Ji a/Xq — 1 
The vortex sheet strength 7 in the previous expression is determined using the Kutta con- 
dition which enforces the physically reasonable condition that the horizontal component of 
the velocity does not diverge at the trailing edge 5 : 

d x (0 7 + (fi nc ) \ X=L = finite (7) 

Substituting (0 EJ) into (J7J) yields the relation 



2ir J 1 \ x — 1 
Multiplying and dividing (jOJ) by the two sides of (JHJ) we obtain 



Xi ' 1 7 dx = Y t + UY x (8) 



P 7 = - WC ^^ C)) PU (Y t + UY X ) (9) 



where 



CM = JT ^7<W A 00 v^T^o (10) 

is the Theodorsen functional j^J which quantifies the unsteadiness of the flow. For example, 
for an airfoil at rest which starts to move suddenly at velocity U, 7 = 6(xq — jUt) corre- 
sponding to the generation of lift due to a vortex that is shed and advected with the fluid. 
Then C = (1 + 2ft/) -1 and we see that asUt/L — > 00, C —>■ 1 , which limit corresponds to the 
realization of the Kutta condition for steady flow |l2|. Adding up the contributions to the 



4 This implies a neglect of any acceleration phase of the vorticity, a reasonable assumption at high Re. 

5 This is tantamount to the statement that that the inclusion of viscosity, no matter how small, will 
regularize the flow in the vicinity of the trailing edge. 
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pressure jump across the filament from the circulatory and non-circulatory flows, we have 
AP = P nc + P 7 , i.e. 

AP = -pUC[y]f (~) (Y t + UY X ) - Lpn (£) Y tt . (11) 
where the dimensionless functions n(s) and f(s) are 



f(s) = 2^— , (12) 
n(s) = 2y/(l - s)s. (13) 

Substituting (fTTj) in (JTJ gives us a single equation of motion for the hydrodynamically driven 
filament 

mY tt = —BY XXXX 

-pUC[i\f{l)(Y t + UY a ) (14) 
-Lpn (§) Y tt . 

with Cpy] determined by (fTUJl . We note that (fT^j) accounts for the unsteady flow past a 
filament of finite length unlike previous studies [8|], and thus includes the effects of vortex 
shedding and fluid added-mass. To make dHj) dimensionless, we scale all lengths with the 
length L of the flag, so that x = sL, Y = rjL, and scale time with the bending time L/Ub, 
where Ug = jj-y/^ is the velocity of bending waves of wavelength 2irL. Then ()14j) may be 
written as 

Mr] TT = -r] ssss - ji8C[j]f(s) (r] T + 5r] s ) (15) 

Here M. = 1 + jj,n(s) where = ~~ = characterizes the added mass effect and the 

parameter 5 = is the ratio of the fluid velocity to the bending wave velocity in the filament. 
We can use symmetry arguments to justify the aerodynamic pressure C[y]f(s) (r) T + dr) s ): 
the term n s arises because the moving fluid breaks the s — > — s symmetry, while the term t] T 
arises because the filament exchanges momentum with the fluid, so that the time reversibility 
t — » —t symmetry is also broken. These two leading terms in the pressure, which could have 
been written down on grounds of symmetry, correspond to a lift force proportional to n s , and 
a frictional damping proportional to r] T . By considering the detailed physical mechanisms, 
we find that the actual form of these terms is more complicated due to the inhomogeneous 
dimensionless functions f(s),n(s). Thus, understanding the flapping instability reduces to 
a stability analysis of the trivial solution rj = of the system (|15I2|) and the determination 
of a transition curve as a function of the problem parameters p, 5. 



Since the free vortex sheet is advected with the flow, the vorticity distribution may 
be written as 7 = j(^-(t — £1) — x ), with (1 + x )/2 denoting the center of vorticity, 
t\ being the time at which shedding occurs; in dimensionless terms reads 7 = 7(25(r — 
t~i) — Xq). Accounting for the oscillatory nature of the flapping instability with an unknown 
frequency uj suggests that an equivalent description of the vorticity distribution is given by 
7 = Ae t( < U)<yT ~ Tl ' > ~ qXoS} where q = to/25 is a non dimensional wave number of the vortex sheet. 
Using the above traveling wave form of the vorticity distribution in (fTU|) we get an expression 
for the Theodorsen function 

c w = c ^ = mntw (16) 

where Hi are Hankel functions of ith order. Substituting the separable form rj(s, r) = £(s)e r7T 
into (|T3|l we get: 

a 2 Mi = -i ssss - C[>y]ti5f(s) K + Ks) (17) 

At the onset of the oscillatory instability, Re(a) = 0, so that a = iu and C[y] is given by 
Then f)17l Ejl constitutes a nonlinear eigenvalue problem for aj,£ given the nonlinear 
uj dependence of the Theodorsen function C(q) = C(u/2S) in We solve the resulting 

system numerically with the AUTO package using a continuation scheme in uj starting 
with a guess for the Theodorsen function C{u/28) = C(0) = 1 . As we shall see later, this 
limit corresponds to the quasi-steady approximation |3|. In Fig. El we show the calculated 
transition curve; when 5 > S c (fi), Re(cr) > with Im(cr) 7^ 0, i.e. an oscillatory instability 
leading to flutter arises. We see that for sufficiently large 5 the filament is always unstable, 
i.e. large enough fluid velocities will always destabilize the elastic filament. As /x ^> 1, the 
added mass effect becomes relatively more important and it is easier for the higher modes 
of the filament to be excited. In Fig. El we show the mode shapes when [i < 1 and \i 3> 1; as 
expected the most unstable mode for /i ^> 1 is not the fundamental mode of the filament. 
We also see that the normalized amplitude of the unstable modes is maximal at the trailing 
edge; this is a consequence of the inhomogeneous functions f(s),n(s) in (|T5j) as well as the 
clamped leading edge and a free trailing edge. 

To further understand the instability, we now turn to a simpler case using the quasi-steady 
approximation [3]. This supposes that the lift forces are slaved adiabatically to those on a 
stationary airfoil with the given instantaneous velocity Y t + UY X , so that C — 1. By assuming 



that the Kutta condition is satisfied instantaneously, we over-estimate the lift forces and thus 
expect to get a threshold for stability that is slightly lower than if C ^ 1. To characterize 
the instability in this situation, we substitute an inhomogeneous perturbation of the form 
i](s,t) = C,(s)e crT into (|15l!2j) and solve the resulting eigenvalue problem to determine the 
growth rate a = a (8, //). In Fig. 2, we show the stability boundary corresponding to the 
quasi-steady approximation. We note that the stability boundary when C ^ 1 accounting 
for vortex shedding corresponds to a higher value of the scaled fluid velocity 5 than that 
obtained using the quasi-steady approximation (7=1, and is a consequence of the quasi- 
steady approximation which over-estimates the lift forces. 
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FIG. 2: Stability diagram for (|15I2I) as a function of the added mass parameter fx and the scaled 
flow velocity 5. The solid line represents the transition curve when vortex shedding is taken into 
account, i.e. C / 1. The dashed line represents the transition curve using the quasi-steady 
approximation where (7=1. In the inset we show the dimensionless wavenumber of the instability 
1 = IS as a function of \i. When g tends to be zero and (7(g) — * 1. The dots correspond to 

experimental data characterizing the transition to flutter in three-dimensional flows past flexible 
sheets of paper the large error bars are a consequence of the variations due to three-dimensional 
effects as well as regions of bistability where both the flapping and stationary state are stable. 

When corresponding to either a fluid of very low density or a filament of very 

high density, Fig. El shows that the corresponding instability occurs for high fluid velocities 
U ^> Ub- Then q —>■ 0, as confirmed in the inset to Fig. |21 Therefore (7(g) = (7(0) = 1 so 
that in this limit the quasi-steady hypothesis is a good approximation. In the limit \i — *■ 0, 
we must have \i8 2 = const so that the aerodynamic pressure which drives the instability 
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FIG. 3: Snapshots of the deflection of the filament rj at the instability threshold for (a) [i = 0.2, 5 
66 and (b) y, = 25, 5 « 6.6. 
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FIG. 4: Spectrum Im(a), Re(a) of the Hamiltonian system ()18I2|) when \x <C 1 (with 5 C 



10.08 > 



The eigenvalues with the smallest absolute value are plotted for 5 = 0.95 c (triangle) and and 
for 5 = 1.15 c (square). We see that instability occurs via a collision and splitting of two pairs of 
eigenvalues along the imaginary axis (indicated by the arrows), and is a signature of a 1:1 resonance 
or a Hamiltonian Hopf bifurcation. 

remains finite. Then the system ()15|) becomes Hamiltonian 6 and may be written as: 



This is because the term breaking time reversal symmetry [i5rj T becomes negligibly small. 
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r] TT = -Vssss - p5 2 f(s)r] s (18) 

The two terms on the right hand side of (|18J) correspond to the existence of two different 
modes of oscillation: (i) that of a flexible filament bending with a frequency that is depen- 
dent on the wavenumber and (ii) that of a rigid filament in the presence of flow-aligning 
aerodynamic forces. In this limiting case, we can clearly see the physical mechanisms at 
work in determining the stability or instability of the filament: small filaments are very stiff 
in bending, but as the filament length becomes large enough for the fluid pressure to excite 
a resonant bending instability the filament starts to flutter. Equivalently, the instability is 
observed when the bending oscillation frequency become of the order of the frequency of 
oscillations of a hinged rigid plate immersed in a flow. To see this quantitatively, we look 
for solutions to ()18I2J) of the form r](s,r) = £(s)e CTT and compute the associated spectrum 
cx(<5). In Fig. HI we show that for 5 < 5 C = 10.08/^//! with < 1, the spectrum lies on 
the imaginary axis as expected, and as 5 > 5 C , the four eigenvalues with smallest absolute 
value collide and split, leading to an instability via a Hamiltonian Hopf Bifurcation or a 1:1 
resonance 3] • 

As fj, ~ 0(1), the effective damping term fj,SCf(s)rj T becomes important, so that the 
spectrum is shifted to the left, i.e. Re(a) < 0. In this case, although the instability is not 
directly related to a resonance, the physical mechanism remains the same, i.e. a competition 
between the destabilizing influence of part of the fluid inertia and the stabilizing influence 
of elastic bending, subject to an effective damping due to fluid motion. This simple picture 
allows us to estimate the criterion for instability by balancing the bending forces -|f with 
the aerodynamic forces IpU 2 ^ so that for a given flow field the critical length of the filament 
above which it will flutter is 

(19) 



IpU 2 

which in dimensionless terms corresponds to 8 ~ 1/ p}^ 2 . Then the typical flapping frequency 
uj is given by balancing filament inertia mu 2 ^ with the aerodynamic forces IpU 2 ^ and leads 
to 

»~M. (20, 



Using typical experimental parameters values from experiments p], we find that L c 
with a frequency uj/2tt = 89 Hz in qualitative agreement with the experimentally observed 
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values L c = 4 cm and a;/27r ~ 50 Hz. In Fig. El we also show the experimental transition 
curve obtained from a recent study on the onset of flutter in paper sheets [5j . The large error 
bars in the experimental data are due to the fact that there is a region of bistability wherein 
both the straight and the flapping sheet are stable. Our linearized theory cannot capture this 
bistability without accounting for the various possible nonlinearities in the system arising 
from geometry. But even without accounting for these nonlinearities, there is a systematic 
discrepancy between our theory and the data which consistently show a higher value of 5 
for the onset of the instability. While there are a number of possible reasons for this, we 
believe that there are two likely candidates: the role of three-dimensional effects and the 
effect of the tension in the filament induced by the Blasius boundary layer, both of which 
would tend to stabilize the sheet and thus push the onset to higher values of S. 

Nevertheless our hierarchy of models starting with the relatively simple Hamiltonian 
picture to the more sophisticated quasi-steady and unsteady ones have allowed us to dissect 
the physical mechanisms associated with flapping in a filament with a finite length and finite 
bending stiffness and account for the added- mass effect, the unsteady lift forces and vortex 
shedding. They also provide a relatively simple criteria for the onset of the instability in 
terms of the scaling laws (|19 [ I20 p . Work currently in progress includes a detailed comparison 
with a two-dimensional numerical simulation and will be reported elsewhere . 
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